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Abstract 

The Gutzwiller wave function for the three-band Hubbard model on a two di- 
mensional Cu02 plane is studied by using the variational Monte Carlo (VMC) 
method in the limit Ud ^ oo, where Ud is the Coulomb repulsion between 
holes on a Cu site. The VMC results for the energy and the fraction of 
d-holes are compared with those of the Rice-Ueda type Gutzwiller approxi- 
mation (GA). The difference between the VMC and the GA results are most 
pronounced at half-filling, and away from half-filling the two results agree well 
in both the hole-doped and the electron-doped cases. The doping dependence 

of the momentum distribution function is also studied. 
PACS numbers: 74.20, 71.30.+h, 71.10.-hx 
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I. INTRODUCTION 



The Gutzwiller wave function (GF) has been used for microscopic investigation of 
strongly correlated electron systems such as the single-band Hubbard model,EHll the pe- 
riodic Anderson model, and the Kondo lattice model. In these last years, the three- 
band Hubbard model on a two dimensional Cu02 plane0 have been studied extensively as 
a model for the high-Tc superconductors. The GF has been also applied to this model with 
the Rice- Ueda0 type Gutzwiller approximation (GA),&@ and with the variational Monte 
Carlo (VMC) method.BHHZl Although the GF is a simple wave function, analytic evaluations 
of expectation values with the GF have been succeeded only in the single-band Hubbard 
model on a one dimensional chain.ill In the VMC method, the expectation values are calcu- 
lated numerically taking the square of the wave function in the first-quantization formalism 
as a probability weight. The merit of the VMC method is that effects of the Gutzwiller 
projection operator which controls the number of doubly occupied sites and spatial correla- 
tions inherent in a starting Fermi sea can be taken into account exactly. These are treated 
with a mean-field like approximation in the GA. Using the VMC method, Yokoyama and 
Shibai have shown that the Brinkman-Rice (BR) transitionil does not occur in the single- 
band Hubbard model. This result has been proved analytically in one dimensional case by 
Metzner and Vollhardt.i 

In previous papers, we have studied the normal state properties in the GF for the three- 
band Hubbard model in the half-filled case,illil where the number of holes per unit cell is 
Uh = 1.0. We have confirmed that the BR transition, which occurs in the GA, is absent 
within the GF also in the three-band Hubbard model by using the VMC method. The 
purpose of this paper is to study effects of doping on the normal state properties. For 
this end, we calculate the energy and the fraction of (i-holes away from half-filling both in 
the hole-dopped {uh > 1.0) and the electron-doped {uh < 1.0) cases using a system which 
contains 6x6 Cu sites. For a comparison, we examine these quantities also by the GA. 
The results show that discrepancies of the VMC and the GA results are most pronounced 
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at half- filling, and away from half- filling the two results agree well when the doping rate \S\, 
which is defined hj = 1.0 + 6, increases. Moreover, we study the momentum distribution 
function n{k) using a 12 x 12 system. Outside the Fermi surface n{k) is nonzero due to the 
strong Coulomb repulsion, and decreases when the wave number k approaches to the jump 
at the Fermi surface. This behavior seems to be one of the characteristics of the GF, and 
was seen also in the GF for the single-band Hubbard modehi'i 

In Sec. H, the GF for the three-band Hubbard model is introduced, and the outline of the 
VMC calculation is explained. The VMC results for the energy and the fraction of (i-holes 
are compared with those of the GA for several doping rates in Sec. |T|. The momentum 
distribution function is studied in Sec. |^ A summary is given in Sec. 0. In the appendix, 



details of the GA calculation are given, and its relation to the auxiliary-boson mean-field 
theory is discussed. 

II. GUTZWILLER WAVE FUNCTION 

The Hamiltonian of the three-band Hubbard model on a two dimensional Cu02 plane is 
given, in a hole notation, by 

n = npd + Ao (Np - Nd) + f/d E « ' (1) 



'Hpd = tpd E ['^Icr ( Pi+S:/2,a + Pi+y/2,a — Pi-x/2,a — Pi-y/2,a ) + h.C. , (2) 
^P = J2 {pl+x/2,aPi+x/2,a + pl+y/2,aPi+y/2,<^) ' (3) 

Nd = E<^ (4) 

i.cr 

where dl^ creates a hole with spin a in the Cu-3 dx2^y2 orbital at site i, nf^ = dj^dia, 
Pi±x/2 cr ipl±y/2 a) creatcs a hole in the 0-px (O-py) orbital, and x {y) is a vector in the x {y) 
direction whose length is equal to the distance between the nearest-neighboring Cu sites. 
The charge transfer energy Aq is written in terms of the on-site energies for a hole on an 
O site Bp and a Cu site Sd as Aq = (^p — £^d)/2 with Ep > Ed- For convenience, the zero 



of energy has been chosen at {Sp + ed)/2 = 0, and we will use the nearest neighboring p-d 
hopping matrix element tp^ as a unit of the energy, i.e., tp^ = 1. In this paper we concentrate 
on the case Ud oo, where Ud is the Coulomb repulsion between holes on a Cu site. In this 
limit doubly occupied states are prohibited on Cu sites. 

The Gutzwiller wave function for the three-band Hubbard model is given by&S 

|^(A)) =Pd|$(A)), (5) 
Pd = n[l-«] > (6) 

i 

where A is the variational parameter, and [^(A)) is a ground state of a noninteracting 
Hamiltonian which contains A as a charge transfer energy 

iC = npd + A (Np - Nd) . (7) 

The wave function |\I'(A)) is an extension of the GF for the single-band Hubbard model0 
and the periodic Anderson model,0'0 and is constructed by projecting out the states with 
doubly occupied Cu sites from the starting Fermi sea (FS), i.e., |<I>(A)). The optimal value of 
A is determined by minimizing the energy for given Hamiltonian parameter Aq. If a Hartree- 
Fock approximation is applied to the Hamiltonian Eq. (||), the on-site energy on a Cu site is 
renormalized due to the Coulomb repulsion, and the charge transfer energy becomes smaller 
than the bare value Aq. This Hartree-Fock picture gives a simple physical interpretation of 
the variational parameter A. The one-particle eigenstates of /C consist of the three bands: 
bonding-, antibonding-, and nonbonding-band. When the hole concentration is < rih < 2, 
the bonding-band, which is the lowest band, is partially filled by holes, and the FS is written 
in the form 

l$(A)) = n 7L|0) , (8) 
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where Vr,k = 2tpd sin(A;^/2), Vyk = 2tpd sm{ky/2), 14 = ^Vj, + V^,,, Sk = \J + V^, 
and is a set of wave numbers in the FS. The operator creates a hole in the bonding- 
band, and dl^, pll, and p^l. are, respectively, the Fourier transforms 

4. = ^^E4e^^«^ (10) 

V^^unit i 
r,^^ - ^ iq{R,+x/2) 

Pka — fJTf 2^Pi+x/2,a^ ' \^^) 

V jVunit i 

r.y\ - ^ \^rit iq{R,+y/2) (-,r,\ 

Pka — nrf 2^Pi+y/2,a^ ' K^"^) 

V^^unit i 

where A^unit is the number of the unit cells each of which contains the three sites Cu02. The 
summation with respect to i (Ri) runs over the Cu sites. 
The energy expectation value, per unit cell , is written as 

E{A) = (H) / N^^,, 

= EpdiA) + Ao [np{A) - nd{A)] , (13) 

where Epd{A) = {Hpd) / N^mt, np{A) = {Np) / N^^^,, nd{A) = {Nd) / N^nit, and 

(O) = <^<^Fl¥» . (14) 
{<I'(A)|*(A)> 

Since the total number of holes Nh is a conserved quantity, np{A) is written in terms of 
nd{A): np{A) = — nd{A) with = / Ny^mt- Thus, the total energy for various Aq 
can be obtained from the data of Epd{A) and nd{A). We calculate the expectation values, 
(O), numerically by using the Monte Carlo (MC) method taking the square of the Slater 
determinant probability weighting. In this method, the constraint of no doubly 

occupied Cu sites and spatial correlations inherent in |$(A)) are treated exactly. Typically, 
20 000-70 000 samples are used for the averages, and each sample is taken every 4-7 MC 
steps after 10 000 MC steps have been discarded for the relaxation, where one MC step 
means a set of NhSNh trials for new configurations. When the acceptance ratio becomes 
small, that is the case for large A, we increase the number of the samples and the sampling 
intervals for keeping the statistical independence of the samples. To estimate the statistical 



fluctuation in the MC calculation, we divide the total samples into 10 groups and, if needed, 
show the maximum and the minimum among them as an error bar. We use a system with 
6x6 unit cells (108 sites) for calculating Epd{A) and nd{A), and use a 12 x 12 system 
for the momentum distribution function. We choose the periodic boundary condition along 
the x-direction, and the antiperiodic boundary condition along the y-direction. With this 
boundary condition, |$(A)) can be determined uniquely in the half-filled case.i 

III. VMC AND GA RESULTS 

In this section, we show the results for the doping dependence of the energy and the 
fraction of cZ-holes. The total number of holes used in the calculation is Nh = 32, 36, 40, 48, 
and 56, and the system size is A^unit = 36. Here the case of rih = 36/36 corresponds to 
half-filling, and the other cases, < 1.0, and, > 1.0, correspond to the electron-doped 
and the hole-doped cases, respectively. We use |tpd| as units of energy. 

In Figs. |I] and |^, the expectation values of the ci-hole fraction ^^(A) and the p-d hopping 
energy Epd{A) are shown as functions of the variational parameter A. The solid circles (•) 
denote the VMC results, and the solid lines denote the corresponding expectation values in 
the FS, i.e., ($(A)|]V,|$(A))/iV,,it or (<l>(A)|7^pd|$(A)) / iV^^if Thus, the differences 
between the circles and the solid lines are due to effects of the projection operator defined 
in Eq. (P). The statistical fluctuations due to the MC calculation are smaller than the size of 
the circles. Both of the results, ^^(A) and Epd{A), show that the difference of the GF from 
the FS becomes large with increasing Uh- Because the doubly occupied states on Cu sites 
are projected out, the fraction of (i-holes is smaller than that in the FS, and nd(A) < 1.0. 
For large A, nd{A) tends to Uh in the electron- doped case (a), and tends to 1.0 in the 
half-filled (b) and the hole-doped cases (c)-(e). The projection operator Pd causes the loss 
of the hopping energy, so that Epd{A) is larger than that in the FS. For large A, Epd{A) 
tends to zero in proportion to l/A.B The hopping energy in the FS becomes minimum at 
A = 0.0, which becomes deep when the hole concentration Uh increases. Correspondingly, 
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Epd{A) becomes minimum at small negative A. In Fig. the value of A which corresponds 
to the minimum of i?pd(A) is plotted as a function of rih- Since the restriction of the p-d 
hopping due to the projection becomes relatively less important when the fraction of 
(i-holes is decreased, the minimum of Epd{A) shifts to the negative value of A. The shift 
becomes large with increasing n^- 

The expectation value of the total energy E{A) is calculated substituting Epd{A), nd{A), 
and np{A) = Uh — nd{A) into Eq. (0). In Fig. ^, E{A) in the half-filled case is shown 
as a function of the variational parameter A for some given Aq. When Aq increases, the 
minimum of E{A) becomes broad, and the value of A which corresponds to the minimum 
becomes large. This is due to the facts that the hopping energy Epd{A) has a minimum at 
A ~ 0, and the potential energy, which is the second term in Eq. (|13]) Aq [np{A) — ^^(A)], 
is a decreasing function of A. The contribution of the potential energy to E{A) in the 
half-filled case is larger than that in the doped cases, which can be seen from the value of 
np{A) — nd{A) in the limit A — ^ oo 



Jim 

A— >oo 



np(A) - HdiA) 



—rih for n/i < 1.0 
72^ — 2 for Ufi > 1.0 

Thus, although the qualitative feature of Fig. H is not changed by doping, the broadening of 
the energy minimum with increasing Aq is most pronounced in the half-filled case. We have 
determined the minimum of E{A) using a least-squares fit to a quadratic function around 
the stationary point. The energy can be determined with reasonable accuracy because an 
error for the energy is proportional to e^, where e is the error for the variational parameter 
A which corresponds to the energy minimum. The error e becomes large when the minimum 
is broad. We will show the error bar due to finite e if it is necessary. In Fig. ^, the optimal 
value of the variational parameter, Aopt, is shown ( function of Aq, and (b) as a 

function of n/j. For small Aq, Aopt is negative and a decreasing function of Uh because 
the total energy is mainly determined by the hopping energy Epd{A). For Aq ^ 2.0, Aopt 
becomes maximum at Uh = 1.0 (b), and the peak becomes sharp with increasing Aq. This is 
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because the total energy is dominated by the potential energy for large Aq. When Aq ^ 1.0, 
the optimal value Aopt is proportional to Ao.H 

In Figs. ^ and |^ the total energy and the d-hole fraction are shown as functions of Aq, 
where the circles (•) denote the VMC results, and the solid lines denote the GA results. 
Details of the GA calculation are given in the appendix. The difference between the VMC 
and the GA energy is large in the half-filled case (b), and becomes small away from half- 
filling both in the electron-doped case (a) and the hole-doped cases (c)-(e). Since the zero of 
energy has been chosen at {ep + ed)/2 = 0, the energy becomes maximum when the fractions 
of p-holes and the d-holes turn out to be the same, i.e., Ud = rip { = nh/2 ). In Fig. ^ the 
error bars for are mainly caused by finite e. In the half-filled case (b), the BR transition 
occurs in the GA at Aq ^ 3.35, and rid becomes unity for Aq > Ag. Although the BR 
transition does not occur within the GF,0 the VMC results for Ud is close to unity for large 
Aq, which means that holes tend to localize at Cu sites. Away from half-filling, the d-ho\e 
fractions with the VMC and the GA agree well. In Fig. the VMC results are plotted as 
functions of Uh for several Aq, where the lines are guide to the eye. For small Aq, the energy 
decreases with increasing n^. However, for Aq ^ 3.0, the VMC energy becomes minimum 
at Uh = 1.0 and seems to show a cusp like behavior. If there is a jump in the slope, the 
system exhibits a phase transition to an insulating state with Ud < 1.0. The d-hole fraction, 
correspondingly, has a peak at nh = 1.0 and becomes close to unity when Aq is large, while 
Ud is an increasing function of Uh for Aq ^ 1.0. We note that beyond the GF the metal- 
insulator transition may occur in the three-band Hubbard model, which has been studied, 
for instance, by Dopf et al with the quantum Monte Carlo method.il 

IV. MOMENTUM DISTRIBUTION FUNCTION 

The momentum distribution function for holes in the bonding-band is defined by 

n{k) = ]-Y.^lllu.). (15) 
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The operator 7]!^ defined in Eq. (^) is written in a linear combination of d and p-lioles, and 
tfie coefficient for eacli orbit depends on tlie variational parameter A. In this definition the 
momentum distribution function in the FS [^(A)) is given by 

^ , ^ I 1.0 if A; G 

($(A)|7L7..|$(A)) = 

[0.0 if A; ^ 

The VMC calculation for n{k) is performed in the real space using the expression which is 
obtained by substituting Eq. and the Fourier transforms Eqs. (|TD|)-(|r^ into Eq. (|T3p. In 
Fig. the momentum distribution function in the half- filled case (a), the hole-doped case 
with rih — 1.0 + 0.833 (b), and the electron-doped case with — 1.0 — 0.833 (c), are plotted 
for several A along the line shown in (d). Because of the antiperiodic boundary condition in 
the ^/-direction, the line is shifted from the usual symmetry directions in the Brillouin zone. 
In the Figure, the value of A is given. The corresponding Hamiltonian parameter Aq is 
estimated from the data obtained in 6 x 6 systems, and is shown in the figure caption. The 
statistical fluctuations in n{k) are smaller than the symbols in the figures. Inside the Fermi 
surface, n{k) is almost flat, and the jump at the Fermi surface becomes small when A (or 
Aq) is large. The jump is finite when 1/A 7^ 0. In the limit A — >■ 00, the jump vanishes 
in the half-filled and the hole-doped cases, while it remains finite in the electron-doped 
case. Nevertheless, in the hole-doped cases, the momentum distribution function for p-holes 
such as {PktPka) ^ jump which remains finite in the limit A — ^ 00. These tendencies 
can be confirmed analytically by expanding •jI^ and |\E'(A)) with respect to 1/A, as it was 
performed in Ref. Outside the Fermi surface, n{k) becomes nonzero due to the Coulomb 
repulsion. The holes excited to the outside are (i-holes because the projection operator Pd 
acts only for ci-holes and commutes with the operators for p-holes. The value of n{k) outside 
the Fermi surface decreases when k approaches to the jump. This behavior seems to be one 
of the characteristics of the Gutzwiller wave function, and was seen also in the GF for the 
single-band Hubbard model.&i 
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V. SUMMARY 



In summary, we have studied the normal-state properties of the three-band Hubbard 
model based on a Gutzwiller wave function, and have compared the VMC results with those 
of the GA. The results for the energy and the fraction of (i-holes show that the discrepancies 
of the VMC and the GA are most pronounced at half-filling. Away from half-filling the 
results with the two methods agree well both in the hole-doped and the electron-doped 
cases. The momentum distribution function for the holes in the bonding-band, n{k), is 
almost flat inside the Fermi surface. Outside the Fermi surface n{k) is finite and decreases 
when k approaches to the jump at the Fermi surface as that in the GF for the single-band 
Hubbard model. 
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APPENDIX A: GUTZWILLER APPROXIMATION 

In this appendix, we demonstrate the calculation of the Rice-Ueda type Gutzwiller ap- 
proximation (GA).0 This approximation was applied to the three-band Hubbard model by 
several authors. Miyake et al and Sarkeiii studied the model with finite Ud, and Mayou 
et a/ii studied the case in which both Ud and Up are oo. Here Up is the Coulomb repulsion 
between holes on an O site. These authors, however, did not investigate the Brinkman-Rice 
transition, which has been studied with the auxiliary-boson method by other authors.0H^ 
The relation between the GA and the auxiliary-boson method is briefly reviewed here. 
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1. Method 



The Gutzwiller wave function in the hmit Ua ^ oo Eq. is generahzed for the model 
with finite Ud by introducing an additional variational parameter g as 

\^{g,A)) = n[l-(l-^7)«]l<f(A)). (Al) 

i 

In the GA, the expectation value for the p-d hopping energy, (^|7ipd|\t')/ (^|^), is replaced 
by that for a noninteracting band with renormalized parameters {^cs\^/(^'Hpd\^cs) ■ Then 
the energy expectation value is given bj"^^ 



EgA = ($eff|/Ceff|$ofr) + UddcuN^nit 

+ (Ao - Aeff) ($eff|(iVp - Nd)\<^es) , (A2) 

where dcn is the fraction of the doubly occupied Cu sites, and l^efr) is a ground state of an 
effective Hamiltonian 

/Ceff =Vqnpd + Aeff {Np - Nd) , (A3) 

with 

ind/2-dcn) ' ' ' 



l-nd + dcn+ ydc 



(A4) 



(n,/2)(l-n,/2) 

In deriving Eq. ([A^), a configurational dependence in the starting Fermi sea |$(A)) defined 
in Eq. (H) is neglected, and a dominant-term approximation is used for dcu and ?7,(i.0'@ The 
condition for the dominant-term with respect to dcu is given by 

^ - {nd/2 - do.? ■ ^^^^ 

Thus the original parameter g is written as a function of Ud and dcu- Since the condition for 
the dominant-term with respect to Ud can not be written in a simple analytic form,0 the 
other original parameter A cannot be determined within the GA. The parameter Aeg in Eqs. 
and ( [A3| ) is introduced in order to treat this condition approximately, and is different 
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from the original parameter A. The value of Aefr and dcu are determined by minimizing the 
energy functional -Ega for given Aq and Ud- 

The effective Hamiltonian /Cefr can be easily diagonalized, and the one-particle eigenstates 
depend only on the ratio of Aeff to the renormalized value of the transfer integral ^/q\tpd\ 

a = Acs/y/q\tpd\ . (A6) 

For < n/j < 2, the d-hole fraction, Ud = {^es\Nd\^c{i) /N^mit, and the ground state energy 
of the effective Hamiltonian, Kefr = {^cslf'^csl^cs) /N^nit, are given by 

1 X- (, ^ Aeff 



and 



1 f'^'^ a 

-rih + / du p{u) , ^ , (A7) 

2 ^-1 Ja^ + 4(1 - u) 



2 



«eff=-^^E VA^g +qV^ 



iVunit 



= -2^q \tpd\ f'^dv p{v)^a^\\{\-v) , (A8) 

where p{y) = J2k^i^ ~ ^'fc)/A''unit, with = {cosk^ + cosky)/2. The parameter is 
determined by the hole concentration as 



nh = 2 du p{u) , (A9) 
and p(i/) is written using the complete elliptic integral of the first kind K{x) as p(z/) = 



2K{\^1 — Substituting Eq. (^) into Eq. (X4), the factor q is written as a function 



of a and dcu- Then other quantities Acg, Kgfj and E'qa are also written as functions of a 
and (icu using Eqs. (|A6|) , ( |A8|) , and Eq. (^) with g = g(Q;,(icu)- When < n/i < 2, the 



conditions for the energy minimum, dEQx/da = and OEqx/ ddQ^ = 0, are obtained as 
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with 



T 



dv p(z/) 



4(1 - y)tla 



(A12) 



2. Results for Ud ^ oo 

In the hmit Ud oo, the doubly occupied states are prohibited on Cu sites, and the 
factor q is written substituting dcn = in to Eq. ( |A4D as 

The optimal value of the ratio a is determined by solving Eq. (|A1(JD. We have calculated 



the integrals in Eqs. (|A7|) - (|A9|) numerically, and obtained the GA results for energy and 
the d-ho\e fraction shown in Figs. ^ and |^. In Fig. |10|, results for q and Aeff are shown 



as functions of Aq for several n^. At half-filling, the Brinkman-Rice transition occurs at 
Aq = Aq 3.35 in units of \tpd\)- The analytic expression for Aq can be obtained by 
expanding Eqj^ with respect to 1/a as 



AS = 2V2^1 + (2/7r)2 \tpd\ . (A14) 

For Aq > Aq, Acff is a two- valued function of Aq when nh = 1.0. The value of Aeff for 
rih = 1 ±0, Afg, can be obtained by substituting Ud = 1.0 and = 0.0 into Eq. ( [A10| ) as 



Ac1r°= AStVA2-A§M/2. (A15) 



For Uh > 1.0, Aeflf tends to zero in the limit Ao/|tpd| oo- This is because the renormal- 
ization factor q tends to zero while the ratio a defined in Eq. (|A6|) has to be smaller than 
a finite value which is determined by Eq. ([A7|) with the condition nd = 1.0. On the other 
hand, in the electron-doped cases there is no upper bound in a, and q > 0.0. Thus, Agff 
increases with Aq for Uh < 1.0. The charge excitation gap egap, defined by a discontinuity 
of the chemical potential, can be also obtained analytically as 
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dE, 



(o) 
GA 



(o) 



GA 



-gap 



dN, 



dNh 



n.=l-0 



(A16) 



2^/A2 



Ac2 
^0 



for Ao < 



for Ao > A^ 



where Eq)^ is the minimized value of -Ega- For Aq > Aq, egap is finite and equal to twice 
of the jump in A^s, i.e., egap = 2 ( A;:^ - A+g° ). 



3. Remarks 

The basic equations presented in this appendix, Eqs. ([A^)- (|A4|) , seems to be equivalent 
to those obtained by Balseiro et al applying a functional integral approach due to Kotliar 
and Ruckensteinil to the three-band Hubbard model. They have, however, neglected the u 
dependence of p(z^), so that their results in the limit Ud oo differ from ours quantitatively. 
The energy functional obtained by Grilli et al within an auxiliary-boson mean-field theory^ 
is also similar to that in the GA. The renormalization factor in their theory, however, is 
different from Eq. ( |A13D , and is equal to the number of the auxiliary-bosons, i.e., q =^ l — Ud- 



Thus, the critical value Aq obtained in auxiliary-boson mean-field theory is a factor -\/2 
smaller than Eq. (|A14| ), which comes from the denominator of Eq. (|A13|) in the limit ^ 1. 
Similar results have also been obtained by Hirashima et al with the exact treatment of 
the local constraint of no doubly occupied Cu sites within the leading order of the 
expansion,ii in which, specifically, the critical value Aq is the same as that of the auxiliary- 
boson mean-field theory. 
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FIGURES 

FIG. 1. The expectation value of the fraction of d holes nd(A) as a function of variational 
parameter A (in units of \tpd\) for several values of the hole concentration Uh- The circles (•) 
denote the VMC results in the Gutzwiller wave function |*(A)), and the solid lines denote the 
d-hole fraction in the Fermi sea |$(A)). 

FIG. 2. The expectation value of the p-d hopping energy, per unit cell, Epd^A) (in units of 
\tpd\) for several rih- The circles (•) denotes the VMC results in the GF, and the lines denote the 
corresponding expectation values in the FS. 

FIG. 3. The value of A which corresponds to the minimum of Ep(i{A) (in units of {tpdD- 

FIG. 4. The expectation value of the total energy, per unit cell, E{A) in the half-filled case 
for Aq = 0.5, 1.0, and 1.5 (in units of \tpd\) with Aq = {sp — £d)/2. The origin of energy is taken 
to be {ep + ed)/2 = 0. 

FIG. 5. The optimal value of the variational parameter Aopt (in units of \tpd\) as a function 
of (a) Ao, and (b) n/j. The lines in (b) are just guide to the eye. 

FIG. 6. Variational results for the total energy, per unit cell, as a function of Aq for several 
Uh (in units of \tpd\)- The circles (•) denote the VMC results, and the solid lines denote the GA 
results. 

FIG. 7. Variational results for the fraction of d-holes rid as a function of Aq (in units of \tpd\)- 
The circles (•) denote the VMC results, and the solid lines denote the GA results. 

FIG. 8. Doping dependence of the VMC results: (a) the total energy per unit cell and (b) the 
cZ-hole fraction, for several values of Aq. The lines are just guide to the eye. 
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FIG. 9. The momentum distribution function for holes in the bonding-band n(fc) defined in 
Eq. (|15|): (a) the half-fihed case {nt = 144/144) with A = 0.1, 1.0, 10.0, and 25.0 (in units of 
\tpd\) which correspond to the optimal value for Aq = 0.3, 1.0, 2.8, and 4.0, respectively, (b) the 
hole-doped case (n/^ = 156/144) with A = 0.1, 1.0, and 10.0 which correspond to the optimal value 
for Aq = 0.4, 1.1, and 4.0, respectively, (c) the electron-doped case (n/i = 132/144) with A = 0.1, 
1.0, and 10.0 which correspond to the optimal value for Aq = 0.2, 0.9, and 3.7, respectively, (d) 
/c-points for a 12 x 12 system with the periodic-antiperiodic boundary condition, where the solid 
and open circles denote, respectively, the occupied and unoccupied states in |$(A)) for rih = 1.0. 

FIG. 10. The GA results for (a) q, and (b) Aggf as functions of Aq (in units of \tpd\) for several 

rih- 
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